bulk-org
Raphael Kubler
12/05/2023
Introduction
Integrated bulk RNA-seq analysis (MSN9, MSN38, WTC11) of cerebral organoid samples in conditions co-cultured (with/without microglia) and stimulation (ctr/LPS/IFNy).
V1.82: added Line-specific analysis
CO vs COiMG
Organoids co-cultured with vs without microglia.
Variance analysis
PCAplot(pca_cor3_mod1, "COiMg", Shape = "Line", PoV.df=PoV_cor3_mod1, pc.1 = 1, pc.2 = 2) #M30 & M33 are outleirs, let's remove them for now
ggplot(pca.df_upgraded,aes_string(x='COiMG',y='Phagocytosis', fill='COiMG')) +
stat_compare_means(method = "t.test", label.y = 10, label.x = 1.3) +
geom_boxplot() +
geom_boxplot(width=0.1, fill="white") +
scale_fill_brewer(palette="Dark2") +
labs(title="Phagocytosis differential PC expression",x="COiMG", y = "PC1 Phagocytosis") +
theme_classic()
Microglia gene heatmap in ctrl samples
#for control only
ctr.meta <- meta[meta$condition %in% 'ctrl',]
ctr.meta$COiMg[ctr.meta$COiMg %in% 'yes'] <- 'COiMg'
ctr.meta$COiMg[ctr.meta$COiMg %in% 'no'] <- 'CO'
#horizontal
ra <- rowAnnotation(
COiMG = ctr.meta$COiMg[ctr.meta$condition %in% 'ctrl'],
col = list(
COiMG = c("COiMg"='tomato','CO'='lightblue')),
show_annotation_name = T,
show_legend = T)
ctr <- t(batch.rem3_mod1[,meta$condition %in% 'ctrl'])
hm_counts <- scale(ctr[,colnames(ctr) %in% gene_panels_subset$Microglia])
ComplexHeatmap::Heatmap(hm_counts,
cluster_rows = T,
cluster_columns = T,
show_row_dend = T,
show_column_dend = T,
show_row_names = F,
show_column_names = F,
left_annotation = ra,
bottom_annotation = HeatmapAnnotation(
text = anno_text(colnames(hm_counts), rot = 300, just = "left",gp=gpar(fontsize=8))),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
name = "Z-score expression")
DEA
Volcano
volcano_plot <- function(res, title = NULL, subtitle = NULL, annotate_by = NULL, type ='ALS', ymax1 = 7.5, ymax2 = 8, xmax1 = -4.5, xmax2 = 4.5){
res <-
mutate(res,
sig = case_when(
padj >= 0.05 ~ "non_sig",
padj < 0.05 & abs(log2FoldChange) < 1 ~ "sig",
padj < 0.05 & abs(log2FoldChange) >= 1 ~ "sig - strong"
)) %>%
mutate(direction = ifelse(log2FoldChange > 0, "up", "down")) %>%
mutate(log2FoldChange = case_when(
log2FoldChange > 10 ~ Inf,
log2FoldChange < -10 ~ -Inf,
TRUE ~ log2FoldChange
)) %>%
mutate(class = paste(sig, direction))
if( type == "ALS"){
xpos <- 0.5
ymax <- ymax1
xlim <- c(xmax1,xmax2)
}else{
xpos <- 0.025
ymax <- ymax2
xlim <- c(-0.042, 0.042)
}
de_tally <- group_by(res, sig, direction, class) %>% tally() %>%
filter(sig != "non_sig") %>%
mutate( position = ifelse(sig == "sig", xpos, 2) ) %>%
mutate(position = ifelse( direction == "down", -1 * position, position)) %>%
mutate(n = formatC(n, format="f", big.mark=",", digits=0))
plot <- res %>%
mutate( pvalue = ifelse( pvalue < 1e-90, Inf, pvalue)) %>% #threshold at 1e16
ggplot(aes(x = log2FoldChange, y = -log10(pvalue))) +
#geom_point(aes(colour = class ), size = 0.5) +
rasterise(geom_point(aes(colour = class ), size = 0.5), dpi = 300) +
scale_colour_manual(values = c("non_sig up" = "gray",
"non_sig down" = "gray",
"sig up" = "#EB7F56",
"sig - strong up" = "#B61927",
"sig down" = "#4F8FC4",
"sig - strong down" = "dodgerblue4"
)) +
theme_bw() +
labs(y = expression(-log[10]~P~value), x = expression(log[2]~"(fold change)"), title = title, subtitle = subtitle) +
guides(colour = FALSE) +
scale_y_continuous(expand = c(0,0), limits = c(0,ymax)) +
theme_bw() +
theme(
plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
panel.border = element_blank(),
axis.ticks = element_line(colour = "black")
) +
geom_text(fontface = "bold", data = de_tally, aes(x = position, y = ymax - 0.5, label = n, colour = class), size = 4) +
scale_x_continuous(limits = xlim)
if(!is.null(annotate_by)){
plot <- plot +
ggrepel::geom_text_repel(
fontface = "italic",
data = filter(res, symbol %in% annotate_by),
aes(x = log2FoldChange, y = -log10(pvalue), label = symbol),
min.segment.length = unit(0, "lines"),
size = 2.3) +
geom_point(
data = filter(res, symbol %in% annotate_by), size = 0.8, colour = "black"
) +
geom_point(aes(colour = class ),
data = filter(res, symbol %in% annotate_by), size = 0.6
)
}
return(plot)
}
results$COiMG$symbol <- rownames(results$COiMG)
volcano_plot(data.frame(results$COiMG), title = 'CO vs COiMg',
annotate_by=c('CD36','CD68','CX3CR1','CSF1R','TREM2','TMEM119','ITGAX'), ymax1 = 30, ymax2 = 31)
Top 50 DEGs heatmap
top50.coimg <- rbind(results$COiMG[order(results$COiMG$log2FoldChange, decreasing = T),][results$COiMG[order(results$COiMG$log2FoldChange, decreasing = T),]$padj < 0.05,][1:25,],
results$COiMG[order(results$COiMG$log2FoldChange, decreasing = F),][results$COiMG[order(results$COiMG$log2FoldChange, decreasing = F),]$padj < 0.05,][1:25,])
top50.names <- rownames(top50.coimg)
ra2 <- rowAnnotation(
COiMG = meta$COiMg,
col = list(
COiMG = c("yes"='tomato','no'='lightblue')),
show_annotation_name = T,
show_legend = T)
ComplexHeatmap::Heatmap(scale(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% top50.names,])),
cluster_rows = T,
cluster_columns = T,
show_row_dend = T,
show_column_dend = T,
show_row_names = F,
show_column_names = F,
left_annotation = ra2,
bottom_annotation = HeatmapAnnotation(
text = anno_text(colnames(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% top50.names,])), rot = 300, just = "left",gp=gpar(fontsize=8))),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
name = "Z-score expression")
Which genes of the top DEGs are cilia-associated (CBC)?
setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)
organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]
coimg.deg <- rbind(results$COiMG[rownames(results$COiMG) %in% na.omit(all_res$COiMG$DEGs[,1]),],
results$COiMG[rownames(results$COiMG) %in% na.omit(all_res$COiMG$DEGs[,2]),])
df <- cbind('DEG'=rownames(coimg.deg),'lFC'=round(coimg.deg$log2FoldChange,2),'Cilia-associated'=rownames(coimg.deg) %in% organoid$CBC)
createDT(data.frame(df), "Status", scrollY=1000)GO
dotplot(all_res$COiMG$GO.up, split="ONTOLOGY", showCategory=10) + facet_grid(ONTOLOGY~., scale="free")
Enrichment
setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)
organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]
#gene set enrichment analysis:
DEG.enrich.res <- matrix(nrow=2*length(colnames(organoid)), ncol=6)
DEG.rep <- NULL
for (i in names(organoid))
DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$COiMG$DEGs[,1]),"Down-regulated"=na.omit(all_res$COiMG$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, organoid[,i], universe=35324))
for (x in colnames(organoid)){
DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(organoid)[,x]), universe=na.omit(rownames(res))))
}
DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
"Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(organoid)
DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
For microglia pathways.
DEG.enrich.res <- matrix(nrow=2*length(colnames(gene_panels_subset)), ncol=6)
DEG.rep <- NULL
for (i in names(gene_panels_subset))
DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$COiMG$DEGs[,1]),"Down-regulated"=na.omit(all_res$COiMG$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, gene_panels_subset[,i], universe=35324))
for (x in colnames(gene_panels_subset)){
DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(gene_panels_subset[,x]), universe=rownames(res)))
}
DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
"Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(gene_panels_subset)
DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
setwd('/Users/kubler01/Documents/R projects/gene lists/ndd')
cgenes <- readxl::read_xlsx('Gene lists for organoid paper.xlsx', col_names = T)
m <- data.frame(t(toupper(na.omit(cgenes$`ADHD GWAS associated genes through FUMA`))[toupper(na.omit(cgenes$`ADHD GWAS associated genes through FUMA`)) %in% rownames(bulkorg_counts2)]), check.names = F)
for(i in names(cgenes)[-1]){
k <- na.omit(cgenes[i])
s <- k[as.matrix(k)[,1] %in% rownames(bulkorg_counts2),]
m <- rbind.fill(data.frame(m),data.frame(t(s)))
}
rownames(m) <- names(cgenes)
m <- t(m)
#gene set enrichment analysis:
ndd.enrich <- matrix(nrow=2*length(colnames(m)), ncol=6)
ndd.rep <- NULL
for (i in colnames(m))
ndd.rep <- c(ndd.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
rownames(ndd.enrich) <- ndd.rep
colnames(ndd.enrich) <- colnames(GSEA.byMod(mod.gl=DEGs, m[,1], universe=35324))
for (x in colnames(m)){
ndd.enrich[gsub("*_up","",gsub("*_down","",rownames(ndd.enrich))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(m)[,x]), universe=na.omit(rownames(res))))
}
ndd.enrich <- data.frame(ndd.enrich)
ndd.enrich$log.q <- -log2(ndd.enrich[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
ndd.enrich.ress <- cbind("Upregulated"=ndd.enrich[c(seq(1, length(rownames(ndd.enrich))-1, by = 2)),][,7],
"Downregulated"=ndd.enrich[c(seq(2, length(rownames(ndd.enrich)), by = 2)),][,7])
rownames(ndd.enrich.ress) <- colnames(m)
print('No developmental DEGs')## [1] "No developmental DEGs"IFNy stimulation
Let’s do this for IFNy stimulation
Variance analysis
PCAplot(pca_cor3_mod1, "condition", Shape = "Line", PoV.df=PoV_cor3_mod1, pc.1 = 1, pc.2 = 2) #M30 & M33 are outleirs, let's remove them for now
ggplot(pca.df_upgraded,aes_string(x='Condition',y='`IFNy`', fill='Condition')) +
#stat_compare_means(method = "t.test", label.y = 20, label.x = 2) +
geom_boxplot() +
geom_boxplot(width=0.1, fill="white") +
scale_fill_brewer(palette="Dark2") +
labs(title="Condition differential PC expression",x="Condition", y = "PC1 IFNy response") +
theme_classic()
DEA
Volcano
results$IFNg$symbol <- rownames(results$IFNg)
volcano_plot(data.frame(results$IFNg), title = 'CTR vs IFNy',
annotate_by=c('CD36','CD68','CX3CR1','CSF1R','TREM2','TMEM119','ITGAX'), ymax1 = 30, ymax2 = 31)
Top 50 DEGs heatmap
top50.ifng <- rbind(results$IFNg[order(results$IFNg$log2FoldChange, decreasing = T),][results$IFNg[order(results$IFNg$log2FoldChange, decreasing = T),]$padj < 0.05,][1:25,],
results$IFNg[order(results$IFNg$log2FoldChange, decreasing = F),][results$IFNg[order(results$IFNg$log2FoldChange, decreasing = F),]$padj < 0.05,][1:25,])
top50.names2 <- rownames(top50.ifng)
ra3 <- rowAnnotation(
Stimulation = meta$condition,
col = list(
Stimulation = c("ctrl"='lightblue','IFNg'='tomato', 'LPS'='darkred')),
show_annotation_name = T,
show_legend = T)
ComplexHeatmap::Heatmap(scale(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% top50.names2,])),
cluster_rows = T,
cluster_columns = T,
show_row_dend = T,
show_column_dend = T,
show_row_names = F,
show_column_names = F,
left_annotation = ra3,
bottom_annotation = HeatmapAnnotation(
text = anno_text(colnames(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% top50.names2,])), rot = 300, just = "left",gp=gpar(fontsize=8))),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
name = "Z-score expression")
Which genes of the top DEGs are cilia-associated (CBC)?
setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)
organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]
IFNg.deg <- rbind(results$IFNg[rownames(results$IFNg) %in% na.omit(all_res$IFNg$DEGs[,1]),],
results$IFNg[rownames(results$IFNg) %in% na.omit(all_res$IFNg$DEGs[,2]),])
df <- cbind('DEG'=rownames(IFNg.deg),'lFC'=round(IFNg.deg$log2FoldChange,2),'Cilia-associated'=rownames(IFNg.deg) %in% organoid$CBC)
createDT(data.frame(df), "Status", scrollY=1000)GO
dotplot(all_res$IFNg$GO.up, split="ONTOLOGY", showCategory=10) + facet_grid(ONTOLOGY~., scale="free")
dotplot(all_res$IFNg$GO.down, split="ONTOLOGY", showCategory=10) + facet_grid(ONTOLOGY~., scale="free")
Enrichment
For organoid cell-types.
setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)
organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]
#gene set enrichment analysis:
DEG.enrich.res <- matrix(nrow=2*length(colnames(organoid)), ncol=6)
DEG.rep <- NULL
for (i in names(organoid))
DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$IFNg$DEGs[,1]),"Down-regulated"=na.omit(all_res$IFNg$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, organoid[,i], universe=35324))
for (x in colnames(organoid)){
DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(organoid)[,x]), universe=na.omit(rownames(res))))
}
DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
"Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(organoid)
DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
For NDD gene lists.
setwd('/Users/kubler01/Documents/R projects/gene lists/ndd')
#gene set enrichment analysis:
ndd.enrich <- matrix(nrow=2*length(colnames(m)), ncol=6)
ndd.rep <- NULL
for (i in colnames(m))
ndd.rep <- c(ndd.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
rownames(ndd.enrich) <- ndd.rep
colnames(ndd.enrich) <- colnames(GSEA.byMod(mod.gl=DEGs, m[,1], universe=35324))
for (x in colnames(m)){
ndd.enrich[gsub("*_up","",gsub("*_down","",rownames(ndd.enrich))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(m)[,x]), universe=na.omit(rownames(res))))
}
ndd.enrich <- data.frame(ndd.enrich)
ndd.enrich$log.q <- -log2(ndd.enrich[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
ndd.enrich.ress <- cbind("Upregulated"=ndd.enrich[c(seq(1, length(rownames(ndd.enrich))-1, by = 2)),][,7],
"Downregulated"=ndd.enrich[c(seq(2, length(rownames(ndd.enrich)), by = 2)),][,7])
rownames(ndd.enrich.ress) <- colnames(m)
ComplexHeatmap::Heatmap(ndd.enrich.ress,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
DEGs$`Down-regulated`[DEGs$`Down-regulated` %in% as.character(na.omit(m[,12]))]## [1] "GRIN2A"For microglia pathways.
DEG.enrich.res <- matrix(nrow=2*length(colnames(gene_panels_subset)), ncol=6)
DEG.rep <- NULL
for (i in names(gene_panels_subset))
DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$IFNg$DEGs[,1]),"Down-regulated"=na.omit(all_res$IFNg$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, gene_panels_subset[,i], universe=35324))
for (x in colnames(gene_panels_subset)){
DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(gene_panels_subset[,x]), universe=rownames(res)))
}
DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
"Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(gene_panels_subset)
DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
Cell-type marker analysis
Cell-types with only 1 gene found back in our expression data carry principal component expression = 0.
Heatmaps of gene & PC expression
setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
cell.types <- as.matrix(readxl::read_xlsx('Typical markers cell types organoids 07182023.xlsx', col_names = F))
PC.covs <- matrix(ncol=length(cell.types[,1]),nrow=length(colnames(batch.rem3_mod1)))
colnames(PC.covs) <- cell.types[,1]
for(i in cell.types[,1]){
pca.cov <- prcomp(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% cell.types[cell.types[,1] %in% i,][-1],]))
PC.covs[,i] <- pca.cov$x[,1]
#if (table(rownames(batch.rem3_mod1) %in% cell.types[cell.types[,1] %in% i,][-1])[[2]] == 1){
# PC.covs[,i] <- t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% cell.types[cell.types[,1] %in% i,][-1],])
#}
}
rownames(PC.covs) <- colnames(batch.rem3_mod1)
ra4 <- rowAnnotation(
Stimulation = meta$condition,
COiMG = meta$COiMg,
col = list(COiMG = c('yes'='black','no'='white'),
Stimulation = c("ctrl"='lightblue','IFNg'='tomato', 'LPS'='darkred')),
show_annotation_name = T,
show_legend = T)
ComplexHeatmap::Heatmap(PC.covs,
cluster_rows = F,
cluster_columns = F,
show_row_dend = T,
show_column_dend = T,
show_row_names = F,
show_column_names = F,
left_annotation = ra4,
bottom_annotation = HeatmapAnnotation(
text = anno_text(colnames(PC.covs), rot = 300, just = "left",gp=gpar(fontsize=8))),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
name = "PC1 expression")
#Let's do that for ctr/ifng and coimg/CO only
#for control only
ctr.meta <- meta[meta$condition %in% 'ctrl',]
#horizontal
ctr.ra <- rowAnnotation(
COiMG = ctr.meta$COiMg,
col = list(
COiMG = c("yes"='tomato','no'='lightblue')),
show_annotation_name = T,
show_legend = T)
ctr.pca <- PC.covs[meta$condition %in% 'ctrl',]
ComplexHeatmap::Heatmap(ctr.pca,
cluster_rows = T,
cluster_columns = T,
show_row_dend = T,
show_column_dend = T,
show_row_names = F,
show_column_names = F,
left_annotation = ctr.ra,
bottom_annotation = HeatmapAnnotation(
text = anno_text(colnames(ctr.pca), rot = 300, just = "left",gp=gpar(fontsize=8))),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
name = "PC1 expression")
ifng.meta <- meta[!meta$condition %in% 'LPS',]
#horizontal
ifng.ra <- rowAnnotation(
COiMG = ifng.meta$COiMg,
Stimulation = ifng.meta$condition,
col = list(COiMG = c('yes'='black','no'='white'),
Stimulation = c("ctrl"='lightblue','IFNg'='tomato', 'LPS'='darkred')),
show_annotation_name = T,
show_legend = T)
ifng.pca <- PC.covs[!meta$condition %in% 'LPS',]
ComplexHeatmap::Heatmap(ifng.pca,
cluster_rows = T,
cluster_columns = T,
show_row_dend = T,
show_column_dend = T,
show_row_names = F,
show_column_names = F,
left_annotation = ifng.ra,
bottom_annotation = HeatmapAnnotation(
text = anno_text(colnames(ctr.pca), rot = 300, just = "left",gp=gpar(fontsize=8))),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
name = "PC1 expression")
#Let's check only choroid plexus & ependymal
cp.ep <- c(na.omit(cell.types[cell.types[,1] %in% c('Chorioid plexus','Ependymal'),][1,][-1]),
na.omit(cell.types[cell.types[,1] %in% c('Chorioid plexus','Ependymal'),][2,][-1]))
ctr.df <- batch.rem3_mod1[,colnames(batch.rem3_mod1) %in% rownames(ctr.meta)]
ComplexHeatmap::Heatmap(scale(t(ctr.df[rownames(ctr.df) %in% cp.ep,])),
cluster_rows = T,
cluster_columns = T,
show_row_dend = T,
show_column_dend = T,
show_row_names = F,
show_column_names = F,
left_annotation = ctr.ra,
bottom_annotation = HeatmapAnnotation(
text = anno_text(colnames(scale(t(ctr.df[rownames(ctr.df) %in% cp.ep,]))), rot = 300, just = "left",gp=gpar(fontsize=8))),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
name = "Z-score expression")
Boxplots
Let’s make boxplots instead
ifng.pca <- PC.covs[!meta$condition %in% 'LPS',]
ifng.pca_long <- reshape::melt(ifng.pca)
colnames(ifng.pca_long) <- c('Sample','Cell.type','PC')
ifng.pca_long$Stimulation <- ifng.meta$condition
ggplot(ifng.pca_long,aes_string(x='Cell.type',y='PC', fill='Stimulation')) +
#stat_compare_means(method = "t.test", label.y = 6) +
geom_boxplot() +
labs(title="Cell-type marker expression",x="Cell-type", y = "PC1 marker expression") +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, vjust = 0.5))
ctr.meta <- meta[meta$condition %in% 'ctrl',]
ctr.pca <- PC.covs[meta$condition %in% 'ctrl',]
coimg.pca_long <- reshape::melt(ctr.pca)
colnames(coimg.pca_long) <- c('Sample','Cell.type','PC')
coimg.pca_long$Coculture <- ctr.meta$COiMg
ggplot(coimg.pca_long,aes_string(x='Cell.type',y='PC', fill='Coculture')) +
#stat_compare_means(method = "t.test", label.y = 6) +
geom_boxplot() +
labs(title="Cell-type marker expression",x="Cell-type", y = "PC1 marker expression") +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, vjust = 0.5))
Dissecting DEA results
CO-IFNy
These are the results of the DEA of CTR vs IFNy in CO (not co-cultured with microglia) only. Cilia-associated genes remain downregulated. We find this as well in CO vs COiMG and more so also in our single-cell data, which did not have any IFNy-stimulated samples. Meaning: this is a result robust to IFNy AND microglia co-culture separately.
df <- cbind('DEG'=rownames(IFNy.res),'lFC'=round(IFNy.res$log2FoldChange,2),'Cilia-associated'=rownames(IFNy.res) %in% organoid$CBC)
createDT(data.frame(df), "Stimulation", scrollY=1000)Volcano
new_IFNy$symbol <- rownames(new_IFNy)
volcano_plot(data.frame(new_IFNy), title = 'ctrl vs IFNy',
annotate_by=c('CD36','CD68','CX3CR1','CSF1R','TREM2','TMEM119','ITGAX'), ymax1 = 30, ymax2 = 31)
CO-COiMg
This is to show that the DEA effect of cilia-associated genes remains present in our bulk-seq analysis when we subset to control samples only and test CO vs COiMg.
df <- cbind('DEG'=rownames(coimg.res),'lFC'=round(coimg.res$log2FoldChange,2),'Cilia-associated'=rownames(coimg.res) %in% organoid$CBC)
createDT(data.frame(df), "COiMG", scrollY=1000)Volcano
new_coimg$symbol <- rownames(new_coimg)
volcano_plot(data.frame(new_coimg), title = 'CO vs COiMg',
annotate_by=c('CD36','CD68','CX3CR1','CSF1R','TREM2','TMEM119','ITGAX'), ymax1 = 30, ymax2 = 31)
Correlation of lFCs
overlap <- rownames(coimg.res[rownames(coimg.res) %in% rownames(IFNy.res),])
cor(coimg.res[rownames(coimg.res) %in% overlap,]$log2FoldChange,IFNy.res[rownames(IFNy.res) %in% overlap,]$log2FoldChange)## [1] 0.718807createDT(cbind('Genes'=overlap,
'COiMg'=coimg.res[rownames(coimg.res) %in% overlap,]$log2FoldChange,
'IFNy'=IFNy.res[rownames(IFNy.res) %in% overlap,]$log2FoldChange), "COiMG", scrollY=1000)By Line analysis
Correlation of DEA results across lines by condition
IFNy
'IFNy_overlap_wtc11Xmsn9:'## [1] "IFNy_overlap_wtc11Xmsn9:"line.res$DEG_IFNy$overlap_wtc11Xmsn9## [1] "CPM" "KIT" "HPGD" "JUP" "KCNJ11" "AKR1C3"
## [7] "LAP3" "MMP25" "CD74" "FAS" "ARAP2" "HERPUD1"
## [13] "PARP12" "OAS1" "FKBP5" "APOL4" "APOL1" "TNFSF13B"
## [19] "GIMAP2" "RIGI" "ACTA2" "LGALS3BP" "DHX58" "CLEC2B"
## [25] "OAS3" "OAS2" "BTN3A3" "CISH" "IFIH1" "GBP3"
## [31] "IFIT3" "IFIT2" "MT2A" "IFI6" "APOBEC3F" "BST2"
## [37] "HELZ2" "GCH1" "TRIM22" "XAF1" "EPSTI1" "PLAAT4"
## [43] "OASL" "SP110" "ALDH1L2" "RTP4" "IL18BP" "IFI44L"
## [49] "IFI44" "PARP9" "HERC6" "PML" "SECTM1" "IFITM3"
## [55] "ITGA9" "NACC2" "SLC43A1" "SLC7A11" "ANKRD22" "MOV10"
## [61] "MX1" "LY6E" "GBP2" "GBP4" "CDCP1" "ERAP2"
## [67] "IFI27" "BATF2" "TAP1" "VAMP5" "MARCHF3" "SAMD9L"
## [73] "PARP10" "CIITA" "DDX60L" "UBA7" "PYCR1" "MX2"
## [79] "USP18" "SOCS1" "IRF7" "IFIT1" "TRIM69" "IFITM1"
## [85] "BTN3A2" "ISG15" "RUFY4" "PLSCR1" "CASP4" "HLA-DOA"
## [91] "HLA-DRA" "HLA-F" "IRF9" "HLA-DPB1" "HLA-DPA1" "HLA-B"
## [97] "CCDC194"print('Correlation of IFNy results between WTC11 & MSN9')## [1] "Correlation of IFNy results between WTC11 & MSN9"cor(line.res$DEG_IFNy[['WTC11']][rownames(line.res$DEG_IFNy[['WTC11']]) %in% line.res$DEG_IFNy$overlap_wtc11Xmsn9,]$log2FoldChange,
line.res$DEG_IFNy[['MSN9']][rownames(line.res$DEG_IFNy[['MSN9']]) %in% line.res$DEG_IFNy$overlap_wtc11Xmsn9,]$log2FoldChange)## [1] 0.2152304'###################'## [1] "###################"'IFNy_overlap_msn9Xmsn38:'## [1] "IFNy_overlap_msn9Xmsn38:"line.res$DEG_IFNy$overlap_msn9Xmsn38## [1] "MMP25" "CD74" "CYP46A1" "TNS1" "CMPK2" "ALDH1L2"
## [7] "IL18BP" "NACC2" "HPGD" "MARCHF3" "CIITA" "RUFY4"
## [13] "HLA-DOA" "HLA-DRA" "HLA-DPB1" "HLA-DPA1" "LAP3" "FAS"
## [19] "ARAP2" "PARP12" "DDX3Y" "SIGLEC1" "OAS1" "APOL4"
## [25] "APOL1" "TNFSF13B" "GIMAP2" "RIGI" "ACTA2" "LGALS3BP"
## [31] "DHX58" "CLEC2B" "OAS3" "OAS2" "BTN3A3" "CISH"
## [37] "IFIH1" "GBP3" "SPP1" "IFIT3" "IFIT2" "MT2A"
## [43] "ID1" "IFI6" "APOBEC3F" "RPS4Y1" "BST2" "HELZ2"
## [49] "GCH1" "TRIM22" "XAF1" "CHI3L1" "EPSTI1" "PLAAT4"
## [55] "OASL" "SP110" "RTP4" "IFI44L" "IFI44" "PARP9"
## [61] "HERC6" "HERC5" "PML" "SECTM1" "IFITM3" "ANKRD22"
## [67] "MOV10" "MX1" "LY6E" "GBP2" "GBP4" "CDCP1"
## [73] "ERAP2" "IFI27" "BATF2" "TAP1" "VAMP5" "SAMD9L"
## [79] "PARP10" "DDX60L" "UBA7" "MX2" "USP18" "SOCS1"
## [85] "IRF7" "IFIT1" "TRIM69" "IFITM1" "BTN3A2" "ISG15"
## [91] "PLSCR1" "CASP4" "HLA-F" "IRF9" "HLA-B" "CCDC194"print('Correlation of IFNy results between MSN38 & MSN9')## [1] "Correlation of IFNy results between MSN38 & MSN9"cor(line.res$DEG_IFNy[['MSN38']][rownames(line.res$DEG_IFNy[['MSN38']]) %in% line.res$DEG_IFNy$overlap_msn9Xmsn38,]$log2FoldChange,
line.res$DEG_IFNy[['MSN9']][rownames(line.res$DEG_IFNy[['MSN9']]) %in% line.res$DEG_IFNy$overlap_msn9Xmsn38,]$log2FoldChange)## [1] 0.1790641'###################'## [1] "###################"'IFNy_overlap_msn9Xmsn38:'## [1] "IFNy_overlap_msn9Xmsn38:"line.res$DEG_IFNy$overlap_msn38Xwtc11## [1] "ITM2A" "HPGD" "GRIN2A" "ERVMER34-1" "CFH"
## [6] "LAP3" "CASP10" "CD38" "NOS2" "MMP25"
## [11] "IL32" "ETV7" "MRC2" "MVP" "NUB1"
## [16] "WAS" "CD74" "SAMD4A" "TYMP" "FAS"
## [21] "CD44" "BTN3A1" "TNFRSF1B" "TNC" "LCP2"
## [26] "ARAP2" "GUCA1A" "EIF2AK2" "PRDM1" "PARP12"
## [31] "OGFR" "LZTS1" "CASP8" "SBNO2" "MTHFD2"
## [36] "SP100" "IFI35" "BCL3" "ASNS" "TBX21"
## [41] "CACNG5" "ARHGAP15" "LAMP3" "SP140" "CEACAM1"
## [46] "IL12RB2" "FYB1" "WDFY1" "CETP" "OAS1"
## [51] "ICAM1" "FLT3LG" "PSME1" "RNF31" "JAK2"
## [56] "HSD3B7" "APOL4" "APOL1" "PCK2" "PSME2"
## [61] "CD40" "TRIB3" "TLDC2" "SAMHD1" "SMCHD1"
## [66] "TNFSF13B" "SLC7A5" "GSDMD" "ZC3HAV1" "CAV2"
## [71] "CAV1" "GIMAP2" "TRIM14" "RIGI" "ACTA2"
## [76] "MAP3K8" "LGALS3BP" "DHX58" "FAM20A" "CPZ"
## [81] "SLC15A3" "CARS1" "CLEC2B" "OAS3" "OAS2"
## [86] "BTN3A3" "NCOA7" "TRIM38" "TENT5A" "ITK"
## [91] "SLC12A7" "BCL6" "CISH" "OTOF" "IFIH1"
## [96] "STAT1" "RNF19B" "FBXO6" "GBP3" "GBP1"
## [101] "IFIT3" "IFIT2" "CD274" "PTK2B" "TRIM25"
## [106] "ABCC11" "TNFSF10" "NRK" "NMI" "BATF3"
## [111] "DNPEP" "ZNFX1" "MT2A" "IRF1" "C3"
## [116] "RBCK1" "IGFLR1" "IFI6" "HELB" "FGL2"
## [121] "STEAP4" "ZFP36" "APOL3" "APOL2" "APOBEC3F"
## [126] "CHAC1" "STARD8" "PLVAP" "BST2" "SAG"
## [131] "HELZ2" "SHFL" "IDO1" "GCH1" "TRIM21"
## [136] "TRIM22" "XAF1" "EPSTI1" "PLAAT4" "DGLUCY"
## [141] "PSRC1" "RSAD2" "IL15RA" "OASL" "PRRG4"
## [146] "PRPH" "SP110" "DOCK10" "ALDH1L2" "PHF11"
## [151] "RTP4" "KLF4" "RNF144B" "IL18BP" "DDX60"
## [156] "CASP1" "SQOR" "IFI44L" "IFI44" "STAMBPL1"
## [161] "PARP9" "HERC6" "CXCL9" "ERP27" "TAPBPL"
## [166] "JDP2" "FBLN5" "WARS1" "PML" "HAPLN3"
## [171] "MARVELD3" "NLRC5" "IRF8" "SECTM1" "PMAIP1"
## [176] "IFITM3" "MOB3C" "FCGR2A" "CYP4V2" "OSMR"
## [181] "TACC1" "NACC2" "SERPING1" "AP1S3" "RASGRP3"
## [186] "ANKRD22" "IFIT5" "LGI4" "GBP5" "L3MBTL4"
## [191] "MOV10" "PIK3AP1" "CLIC2" "ART3" "MAP3K7CL"
## [196] "UBE2L6" "MX1" "TNFRSF14" "MEGF11" "SLAMF8"
## [201] "IFNAR2" "CBR3" "C1R" "ADAR" "LY6E"
## [206] "CXCL16" "PDPN" "GBP2" "GBP4" "VCAM1"
## [211] "ATF3" "CTSS" "FZD5" "IFI16" "CDCP1"
## [216] "DTX3L" "LIPH" "TMEM144" "IL15" "ERAP1"
## [221] "ERAP2" "TLR3" "CASP7" "OTOGL" "IFI27"
## [226] "C2" "PLEKHF1" "TVP23A" "B2M" "TBC1D2B"
## [231] "NOD2" "PRRT2" "BATF2" "IRF2" "FILIP1L"
## [236] "TAP1" "MLKL" "STAT3" "VAMP5" "LGALS9"
## [241] "ATF5" "CXCL10" "CXCL11" "PLEKHA2" "PROKR1"
## [246] "TRIM56" "TCAF2" "STAT2" "TMEM51" "ISG20"
## [251] "CEBPB" "SLFN11" "MYD88" "PARP14" "PARP15"
## [256] "C1QB" "C1QA" "CD7" "RNF213" "MARCHF3"
## [261] "DES" "DDIT3" "A2M" "SAMD9L" "PARP10"
## [266] "EXOC3L1" "CIITA" "DDX60L" "UBA7" "C1S"
## [271] "MX2" "ASCL2" "FAM162B" "CSF1" "SOCS3"
## [276] "USP18" "MAFF" "TNFAIP2" "SOCS1" "SP140L"
## [281] "IRF7" "IFIT1" "TRIM69" "IFITM1" "BTN3A2"
## [286] "ISG15" "EIF4EBP1" "RUFY4" "PLSCR1" "CALHM6"
## [291] "INSYN2A" "ANKRD34B" "HLA-DRB1" "ARID5A" "CASP4"
## [296] "ACSL5" "GSTK1" "MAP1LC3C" "TEAD4" "TGM2"
## [301] "HLA-DOA" "HLA-DMA" "PSMB8" "TAP2" "HLA-DRA"
## [306] "CARD16" "LST1" "HLA-C" "HLA-E" "HLA-F"
## [311] "PSMB10" "SAMD9" "ATP10A" "HLA-A" "UBD"
## [316] "IRF9" "SMTNL1" "IFI30" "CEBPD" "APOL6"
## [321] "HLA-DPB1" "C4B" "PATL2" "HLA-DPA1" "CYP21A2"
## [326] "TAPBP" "HLA-B" "KIAA0040" "OR2I1P" "APOBEC3G"
## [331] "PSMB9" "HLA-DOB" "HLA-DMB" "CFB" "APOBEC3D"
## [336] "APOBEC3C" "C4A" "TMEM158" "PRSS51" "TRIL"
## [341] "ZNF878" "TRIM34" "CCDC194" "CCL5" "PRAG1"print('Correlation of IFNy results between MSN38 & WTC11')## [1] "Correlation of IFNy results between MSN38 & WTC11"cor(line.res$DEG_IFNy[['MSN38']][rownames(line.res$DEG_IFNy[['MSN38']]) %in% line.res$DEG_IFNy$overlap_msn38Xwtc11,]$log2FoldChange,
line.res$DEG_IFNy[['WTC11']][rownames(line.res$DEG_IFNy[['WTC11']]) %in% line.res$DEG_IFNy$overlap_msn38Xwtc11,]$log2FoldChange)## [1] 0.9095827LinexCondition
for (i in c('FOXJ1',
'CLDN5',
'IDO1',
'GBP5',
'MMP10',
'SOX2',
'RBFOX3')){
p <- ggplot(pca.df_upgraded,aes_string(x='LinexCondition',y=i, fill='LinexCondition')) +
#stat_compare_means(method = "t.test", label.y = 2000, label.x = 1.3) +
geom_boxplot() +
geom_boxplot(width=0.1, fill="white") +
scale_fill_brewer(palette="Dark2") +
labs(title=paste0("LinexCondition ", i, " expression"),x="LinexCondition", y = paste0(i,' expression')) +
theme_classic()
plot(p)
}
COiMg
'IFNy_overlap_wtc11Xmsn9:'## [1] "IFNy_overlap_wtc11Xmsn9:"line.res$DEG_COiMg$overlap_wtc11Xmsn9## [1] "IL1RN" "IFI44" "GIPC3" "ADAMTSL2"print('Correlation of IFNy results between WTC11 & MSN9')## [1] "Correlation of IFNy results between WTC11 & MSN9"cor(line.res$DEG_COiMg[['WTC11']][rownames(line.res$DEG_COiMg[['WTC11']]) %in% line.res$DEG_COiMg$overlap_wtc11Xmsn9,]$log2FoldChange,
line.res$DEG_COiMg[['MSN9']][rownames(line.res$DEG_COiMg[['MSN9']]) %in% line.res$DEG_COiMg$overlap_wtc11Xmsn9,]$log2FoldChange)## [1] 0.9463698'###################'## [1] "###################"'IFNy_overlap_msn9Xmsn38:'## [1] "IFNy_overlap_msn9Xmsn38:"line.res$DEG_IFNy$overlap_msn9Xmsn38## [1] "MMP25" "CD74" "CYP46A1" "TNS1" "CMPK2" "ALDH1L2"
## [7] "IL18BP" "NACC2" "HPGD" "MARCHF3" "CIITA" "RUFY4"
## [13] "HLA-DOA" "HLA-DRA" "HLA-DPB1" "HLA-DPA1" "LAP3" "FAS"
## [19] "ARAP2" "PARP12" "DDX3Y" "SIGLEC1" "OAS1" "APOL4"
## [25] "APOL1" "TNFSF13B" "GIMAP2" "RIGI" "ACTA2" "LGALS3BP"
## [31] "DHX58" "CLEC2B" "OAS3" "OAS2" "BTN3A3" "CISH"
## [37] "IFIH1" "GBP3" "SPP1" "IFIT3" "IFIT2" "MT2A"
## [43] "ID1" "IFI6" "APOBEC3F" "RPS4Y1" "BST2" "HELZ2"
## [49] "GCH1" "TRIM22" "XAF1" "CHI3L1" "EPSTI1" "PLAAT4"
## [55] "OASL" "SP110" "RTP4" "IFI44L" "IFI44" "PARP9"
## [61] "HERC6" "HERC5" "PML" "SECTM1" "IFITM3" "ANKRD22"
## [67] "MOV10" "MX1" "LY6E" "GBP2" "GBP4" "CDCP1"
## [73] "ERAP2" "IFI27" "BATF2" "TAP1" "VAMP5" "SAMD9L"
## [79] "PARP10" "DDX60L" "UBA7" "MX2" "USP18" "SOCS1"
## [85] "IRF7" "IFIT1" "TRIM69" "IFITM1" "BTN3A2" "ISG15"
## [91] "PLSCR1" "CASP4" "HLA-F" "IRF9" "HLA-B" "CCDC194"print('Correlation of IFNy results between MSN38 & MSN9')## [1] "Correlation of IFNy results between MSN38 & MSN9"cor(line.res$DEG_COiMg[['MSN38']][rownames(line.res$DEG_COiMg[['MSN38']]) %in% line.res$DEG_COiMg$overlap_msn9Xmsn38,]$log2FoldChange,
line.res$DEG_COiMg[['MSN9']][rownames(line.res$DEG_COiMg[['MSN9']]) %in% line.res$DEG_COiMg$overlap_msn9Xmsn38,]$log2FoldChange)## [1] 0.6876617'###################'## [1] "###################"'IFNy_overlap_msn9Xmsn38:'## [1] "IFNy_overlap_msn9Xmsn38:"line.res$DEG_COiMg$overlap_msn38Xwtc11## [1] "BCAS1" "DBX1" "FGR" "TMEM176A" "ITGAL" "TRAF3IP3"
## [7] "STAB1" "CD4" "BTK" "MRC2" "TYROBP" "ALOX5"
## [13] "SLC11A1" "RUNX3" "SLAMF7" "TNFRSF1B" "VCAN" "MSR1"
## [19] "CAPG" "ADAM28" "LCP2" "COL9A2" "ELN" "TBXAS1"
## [25] "GNA15" "CD84" "SPI1" "APBB1IP" "EDN1" "SP140"
## [31] "FYB1" "LAT2" "NOX4" "ADAMTS2" "SIGLEC1" "RGS1"
## [37] "LYZ" "UNC13D" "TREM2" "IL12RB1" "CYTH4" "HMOX1"
## [43] "NCF4" "CSF2RB" "HCK" "MXRA5" "ACP5" "PYCARD"
## [49] "AQP9" "IL4I1" "RASAL3" "SIGLEC8" "CD33" "PIK3CG"
## [55] "TFEC" "TMEM176B" "ENG" "ABI3" "COL1A1" "IL10RA"
## [61] "SLC15A3" "C11orf21" "BIN2" "ARHGDIB" "PTPN6" "MDFI"
## [67] "LOX" "CD86" "CYTIP" "ITGB6" "PROC" "PLEK"
## [73] "NCF2" "CD48" "PADI2" "SPP1" "CSF3R" "TGFBI"
## [79] "CCRL2" "TMIGD3" "SASH3" "SRGN" "BHLHE41" "ARHGAP9"
## [85] "NCKAP1L" "PMEPA1" "C3" "EVI2A" "ADGRE2" "FGL2"
## [91] "RAC2" "IRF5" "WDFY4" "CD68" "SIGLEC9" "APOC1"
## [97] "LSP1" "COL5A1" "FIBCD1" "CARD6" "ALOX5AP" "CHI3L1"
## [103] "CHIT1" "PRAM1" "GIMAP6" "GIMAP4" "ADAMDEC1" "CD180"
## [109] "PTPN22" "FST" "DOCK2" "SLC37A2" "HAVCR2" "SDS"
## [115] "CD36" "NT5E" "DYSF" "NPL" "LCP1" "IL1RN"
## [121] "TLR4" "SLCO2B1" "CASP1" "ITGA11" "PLCB2" "PARVG"
## [127] "ACVRL1" "GALNT6" "GPR65" "BCL2A1" "GOLGA6D" "ITGAX"
## [133] "PIK3R5" "VAV1" "MYO1F" "CD53" "FCGR2A" "LEFTY2"
## [139] "PTPN7" "OSBPL10" "TM4SF19" "PLA2G7" "DOK2" "CDKN2B"
## [145] "GSN" "RGS10" "ST14" "TAGLN" "FERMT3" "FCGR1A"
## [151] "BMP6" "SAMSN1" "SLC7A7" "VSIG4" "TDRD9" "CD109"
## [157] "SLC34A2" "TNFRSF14" "NCF1" "SLAMF8" "FCER1G" "C1QC"
## [163] "RUNX1" "GOLGA6A" "ITGB2" "IL6R" "LAPTM5" "FCMR"
## [169] "HPGDS" "CTSS" "S100A11" "S100A9" "CCKAR" "FBLN2"
## [175] "CCR1" "TMEM144" "TAGAP" "DCSTAMP" "SYK" "ALDH1A1"
## [181] "CYBB" "HTRA1" "SPINT1" "PLD4" "MFAP4" "C15orf48"
## [187] "MS4A7" "MS4A14" "LAIR1" "LDLRAD4" "CXCL8" "CD52"
## [193] "ITGAM" "FPR1" "GPR34" "C3AR1" "SYNPO" "KLHL6"
## [199] "C1QB" "C1QA" "OLR1" "TLR10" "TLR6" "LRRC25"
## [205] "NUPR1" "MSC" "SSC5D" "HCLS1" "CLDN7" "UBA7"
## [211] "BGN" "CSF1R" "TMEM119" "APOBR" "EVI2B" "CD300LF"
## [217] "FFAR4" "ARHGAP30" "TRPV2" "AMTN" "DPYD" "TLR7"
## [223] "SPN" "ADAMTSL2" "RCSD1" "SUCNR1" "FCGR3A" "LST1"
## [229] "STK38L" "APOC2" "HLA-DMB" "CLEC5A" "CLEC19A" "PECAM1"
## [235] "RAB7B" "CCL3" "ADORA3"print('Correlation of IFNy results between MSN38 & WTC11')## [1] "Correlation of IFNy results between MSN38 & WTC11"cor(line.res$DEG_COiMg[['MSN38']][rownames(line.res$DEG_COiMg[['MSN38']]) %in% line.res$DEG_COiMg$overlap_msn38Xwtc11,]$log2FoldChange,
line.res$DEG_COiMg[['WTC11']][rownames(line.res$DEG_COiMg[['WTC11']]) %in% line.res$DEG_COiMg$overlap_msn38Xwtc11,]$log2FoldChange)## [1] 0.8966225LinexCOiMG
#LinexCOiMG
for (i in c('FOXJ1',
'CLDN5',
'IDO1',
'GBP5',
'MMP10',
'SOX2',
'RBFOX3')){
p <- ggplot(pca.df_upgraded,aes_string(x='LinexCOiMG',y=i, fill='LinexCOiMG')) +
#stat_compare_means(method = "t.test", label.y = 2000, label.x = 1.3) +
geom_boxplot() +
geom_boxplot(width=0.1, fill="white") +
scale_fill_brewer(palette="Dark2") +
labs(title=paste0("LinexCOiMG ", i, " expression"),x="LinexCOiMG", y = paste0(i,' expression')) +
theme_classic()
plot(p)
}
Interaction analysis
Interaction analysis of co-culture*IFNy
What happens if you co-culture AND stimulate with IFNy?
Variance analysis
pca_cor3_mod1$interaction <- as.factor(pca_cor3_mod1$COiMg):as.factor(pca_cor3_mod1$condition)
PCAplot(pca_cor3_mod1, "interaction", Shape = "Line", PoV.df=PoV_cor3_mod1, pc.1 = 1, pc.2 = 2,
colors = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"Spectral")))(6)) #M30 & M33 are outleirs, let's remove them for now
ggplot(pca.df_upgraded,aes_string(x='Condition',y='`IFNy`', fill='Condition')) +
#stat_compare_means(method = "t.test", label.y = 20, label.x = 2) +
geom_boxplot() +
geom_boxplot(width=0.1, fill="white") +
scale_fill_brewer(palette="Dark2") +
labs(title="Condition differential PC expression",x="Condition", y = "PC1 IFNy response") +
theme_classic()
DEA
Volcano
results$COiMGxIFNg$symbol <- rownames(results$COiMGxIFNg)
volcano_plot(data.frame(results$COiMGxIFNg), title = 'CTR vs IFNy',
annotate_by=c('CD36','CD68','CX3CR1','CSF1R','TREM2','TMEM119','ITGAX'), ymax1 = 10, ymax2 = 11)
DEA boxplots
COiMGxCondition
for (i in c('FOXJ1',
'CLDN5',
'IDO1',
'GBP5',
'MMP10',
'SOX2',
'RBFOX3')){
p <- ggplot(pca.df_upgraded,aes_string(x='COiMGxCondition',y=i, fill='COiMGxCondition')) +
#stat_compare_means(method = "t.test", label.y = 2000, label.x = 1.3) +
geom_boxplot() +
geom_boxplot(width=0.1, fill="white") +
scale_fill_brewer(palette="Dark2") +
labs(title=paste0("COiMGxCondition ", i, " expression"),x="COiMGxCondition", y = paste0(i,' expression')) +
theme_classic()
plot(p)
}
